COMMENT ⊗   VALID 00008 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	% This is a sample WEB program written by DEK in March, 1983.
C00003 00003	@* Introduction.
C00007 00004	@* Arithmetic.
C00013 00005	@* Representing the data.
C00017 00006	@* Row operations.
C00023 00007	@* The main program.
C00024 00008	@* Index.
C00025 ENDMK
C⊗;
% This is a sample WEB program written by DEK in March, 1983.
\def\topofcontents{\null
	\vfill
	\centerline{\titlefont A routine for vector processing}
	\vfill}
@* Introduction.
The purpose of this program is to find linear dependencies that exist
between $m$@@given integer vectors of dimension@@$n$, and to determine
the ``null space'' of those vectors. For example, suppose that $n=4$ and
that $m=5$ vectors are given: $a=(1,1,0,-1)$, $b=(1,0,1,0)$, $c=(1,1,0,1)$,
$d=(0,1,-1,0)$, $e=(2,1,1,-1)$. The program will discover that
$d={1\over2}a-b+{1\over2}c$ and that $e=a+b$; and it will also print the
vector $v=(-1,1,1,0)$ which has the property that $v\cdot a=v\cdot b=
v\cdot c=v\cdot d=v\cdot e=0$.

In general, suppose that the given set of vectors has rank@@$r$; then there
will be $m-r$ dependencies, and $n-r$ vectors in the null space.

@ The input data is given on file |input| in the following form:
First comes the value of $n$, on a line by itself;
then comes $m$ lines of the form $c\,a_1\ldots a_n$, where $c$ is a
one-character label for the vector, and where $a_1\ldots a_n$ are the
components; then comes the end of the file.

For example, the input data mentioned above would appear in six lines:
$$\vbox{\halign{\tt#\cr
4\cr
a 1 1 0 -1\cr
b 1 0 1 0\cr
c 1 1 0 1\cr
d 0 1 -1 0\cr
e 2 1 1 -1\cr
}}$$

The |output| file will contain a copy of the input data, together with
any dependencies found, followed by the null space data. It will look like
this in our example:
$$\vbox{\halign{\tt#\cr
a=(1,1,0,-1)\cr
b=(1,0,1,0)\cr
c=(1,1,0,1)\cr
d=(0,1,-1,0)\cr
 =(1/2)a-b+(1/2)c.\cr
e=(2,1,1,-1)\cr
 =a+b.\cr
The following vectors generate the null space:\cr
(-1,1,1,0)\cr
}}$$

@ The program is written entirely in standard {\tt PASCAL}. (Or, at least
the author tried to do so.)

@p program vector(input,output);
const @!max_n=50; {upper bound on $n$}
type @<Types in the outer block@>@;@/
var @<Global variables@>

@ In one place the author resorts to |goto| statements. The following
symbolic labels make this slightly less sinful.

@d exit=30 {go here to return from a procedure}
@d found=40 {go here when you've found what you were looking for}
@* Arithmetic.
All of the arithmetic in this program is done with fractions, since the
numbers will be small in all intended applications, and since this
avoids problems of numerical instability. We assume that the numerator
and denominator will not exceed the range of normal integers.

In each variable of type |fraction|, there are two fields, |num| and |denom|.
The values of these fields will be integers having no common factor,
and |denom| will be positive.

It follows that a fraction is an integer if and only if |denom=1|, and it
is zero if and only if |num=0|.

@d nonzero(#)==(#.num≠0) {test if a fraction is nonzero}
@d nonint(#)==(#.denom>1) {test if a fraction is not a whole number}

@<Types...@>=
@!fraction=record @!num:integer; {the numerator of this fraction}
	@!denom:integer; {the denominator, which is always positive}
	end;

@ For simplicity, all arithmetic is implemented with respect to
an accumulator, which is a global |fraction|.

@d load(#)==acc←#
@d load_negative(#)==begin acc.num←-#.num; acc.denom←#.denom;
	end
@d store(#)==#←acc
@d store_int_end(#)==#;@+end
@d store_int(#)==begin #.denom←1;#.num←store_int_end
	{one says |store_int(v)(x)|}

@<Glob...@>=
@!acc:fraction; {the input and output of arithmetic operations}

@ Since we work with fractions, it should not be surprising that
we are going to use Euclid's algorithm.

@p function gcd(@!u,@!v:integer):integer;
var w:integer; {temporary storage}
begin while v≠0 do
	begin w←u; u←v; v←w mod v;
	end;
gcd←abs(u);
end;

@ The following procedures implement the methods described in Section@@4.5.1
of {\sl The Art of Computer Programming}. We use the abbreviations
|un|, |ud|, |vn|, and@@|vd| for what the book calls $u$, $u'$, $v$, and@@$v'$.

@d un==acc.num
@d ud==acc.denom
@d vn==x.num
@d vd==x.denom

@p procedure add(var x:fraction); {set |acc←acc+x|}
var d1,@!d2,@!t:integer; {temporary storage}
begin d1←gcd(ud,vd);
if d1=1 then
	begin un←un*vd+ud*vn; ud←ud*vd;
	end
else	begin t←un*(vd div d1)+vn*(ud div d1);
	d2←gcd(t,d1);
	un←t div d2; ud←(ud div d1)*(vd div d2);
	end;
end;
@#
procedure multiply_by(var x:fraction); {set |acc←acc*x|}
var d1,@!d2:integer; {temporary storage}
begin d1←gcd(un,vd); d2←gcd(ud,vn);
un←(un div d1)*(vn div d2);
ud←(ud div d2)*(vd div d1);
end;
@#
procedure divide_by(var x:fraction); {set |acc←acc/x|, given |x≠0|}
var d1,@!d2:integer; {temporary storage}
begin d1←gcd(un,vn); d2←gcd(ud,vd);
un←(un div d1)*(vd div d2);
ud←(ud div d2)*(vn div d1);
if ud<0 then
	begin un←-un; ud←-ud;
	end;
end;

@ Here's a procedure that outputs a fraction.

@p procedure write_fraction(var x:fraction);
begin write(x.num:1);
if nonint(x) then write('/',x.denom:1);
end;

@ A here's a similar one that is used to output the coefficient of a vector,
when that coefficient is a nonzero fraction.

@p procedure write_coef; {output |acc| as a coefficient}
begin if acc.num<0 then write('-')
else if initial_plus then write('+');
if nonint(acc) then
	begin write('('); write(abs(acc.num):1); write('/');
	write(acc.denom:1); write(')');
	end
else if abs(acc.num)≠1 then write(abs(acc.num):1);
initial_plus←true;
end;

@ @<Glob...@>=
@!initial_plus:boolean; {should |write_coef| write a plus?}
@* Representing the data.
Now let's consider the main data structures that are used.

@ There's an array |name| of one-character vector labels, copied from the
input. (We remember the name only if the vector was independent of all
that preceded it.)

@<Glob...@>=
@!name:array [0..max_n] of char; {vector labels}

@ Because of the nature of the algorithm that we will develop, it turns
out to be convenient to work with vectors that are ``split at |i|''.
This means that we interpret the entries |v[j]| for |j≤i| differently
from the entries for |j>i|. Such a vector will satisfy the condition
$$(0,\ldots,0,1,v[i+1],\ldots,v[n])=\sum↓{1\le j\le i}v[j]\cdot|name|[j]$$
where |name[j]| stands for the input vector of that name.
In other words, the positions to the left of@@|i| tell us how the
vector was formed; the vector itself will have nothing but zeroes
to the left of position@@|i|, and position@@|i| itself will be@@1.

@<Glob...@>=
@!i,@!j,@!k:integer; {miscellaneous indices}

@ Most of the data appears in an $n\times n$ array |A| of fractions.
Row@@|i| of this array---i.e., the vector $(A[i,1],\ldots,A[i,n])$---
is either an element of the null space or a vector that is split at@@|i|
as described above. An auxiliary boolean array called |null_space|
distinguishes these two cases: If |null_space[i]| is |true|,
row@@|i| is an element of the null space. Furthermore |A[i,i]=1| in this
case, and |A[i,j]=0| for |j>i|; and |A[i,j]≠0| for |j<i| only
if |null_space[j]| is |false|.

@<Glob...@>=
@!n:1..max_n; {the size of the vectors}
@!A:array[1..max_n,1..max_n] of fraction; {the main row vector}
@!null_space:array[1..max_n] of boolean; {the types of rows in |A|}

@ Here's how those data structures can be initialized so that the
stated properties hold: We simply set |A| to the identity matrix,
and mark all of its rows as belonging to the null space.

@<Initialize the data structures@>=
read_ln(n);
for i←1 to n do
	begin null_space[i]←true;
	for j←1 to n do
		if i=j then store_int(A[i,j])(1)
		else store_int(A[i,j])(0);
	end
@* Row operations.
The main point of those data structures is that we can maintain them
when a new vector is read from the input.

The algorithm below reads a new vector into a working array called@@|v|,
which is a special split array: It has an extra component |v[0]|,
which records its dependence on |name[0]|, which is the label of the
newly input vector. At most times during the processing, the relation
$$(0,\ldots,0,v[i],v[i+1],\ldots,v[n])=\sum↓{0\le j<i}v[j]\cdot|name|[j]
  \eqno({*})$$
will be valid, while we will have
$$(0,\ldots,0,1,v[i+1],\ldots,v[n])=\sum↓{0\le j\le i}v[j]\cdot|name|[j]
  \eqno({**})$$
as in our previous definition at the point where |null_space[i]| is tested.
The reason for having two slightly different conventions is, of course,
that transitions occur at component |i| when |i| increases.

@<Glob...@>=
@!v:array[0..max_n] of fraction; {an auxiliary array}

@ Here is the procedure that does most of the work.

@p procedure new_vector; {absorbs a new vector from the input}
label exit,found;
var i,@!j,@!k:integer; {miscellaneous indices and temporary storage}
@!c:fraction; {dot product used to update the null matrix}
begin @<Read a new vector into |v|, and set |i←0|@>;
while i<n do
	begin i←i+1; {now $(*)$ holds}
	if nonzero(v[i]) then
		begin @<Divide |v| by |v[i]|@>;
		store_int(v[i])(0); {now $(**)$ holds}
		if null_space[i] then goto found;
		@<Subtract row |i| from |v|@>;
		end;
	end;
@<Output the dependency that was just discovered@>;
goto exit;
found:@<Update |A| for a new independent vector |v|@>;
exit:end;

@ @<Read a new vector...@>=
begin read(name[0]); write(name[0],'=(');
for i←1 to n do
	begin read(k); store_int(v[i])(k);
	if i>1 then write(',');
	write(k:1);
	end;
read_ln; write_ln(')');
i←0; store_int(v[0])(1);
end

@ @<Divide |v|...@>=
for j←0 to n do if j≠i then
	begin load(v[j]); divide_by(v[i]); store(v[j]);
	end

@ When the following code is executed, both |v| and $(A[i,1],\ldots,A[i,n])$
are split at@@|i|. Hence we can subtract them, obtaining a vector that
satisfies $(**)$, except that the |i|th component on the left-hand side
is@@0 instead of@@1.

@ @<Subtract row |i|...@>=
for j←1 to n do
	begin load_negative(A[i,j]); add(v[j]); store(v[j]);
	end

@ @<Output the dependency...@>=
begin write(' ='); initial_plus←false;
for i←1 to n do if nonzero(v[i]) then
	begin load_negative(v[i]); divide_by(v[0]); write_coef; write(name[i]);
	end;
write_ln('.');
end

@ When a new independent vector has been found, we have |v[i]=0|; so we
can move |v[0]| into position@@|i|, obtaining a row vector that is
split at@@|i| as required in matrix@@|A|. Before we store that new row,
we must update the remaining null space vectors, by subtracting
|c|@@times row@@|i| from row@@|j|, where |c| is the dot product of
row@@|j| with the new vector.

@<Update |A|...@>=
begin for j←i+1 to n do if null_space[j] then
	begin @<Compute the dot product, |c|@>;
	@<Subtract |c| times row |i| from row |j|@>;
	end;
name[i]←name[0]; v[i]←v[0]; null_space[i]←false;
for j←1 to n do A[i,j]←v[j];
end

@ At this point we know that row@@|j| has zero in position@@|i|, and in
positions greater than@@|j|; furthermore the new vector is zero in positions
less than@@|i|. Therefore we need only consider positions |i+1| through@@|j|.

@<Compute the dot product...@>=
begin store_int(c)(0);
for k←i+1 to j do if nonzero(A[j,k]) then
	begin load(A[j,k]); multiply_by(v[k]); add(c); store(c);
	end;
end

@ @<Subtract |c| times...@>=
for k←1 to i do
	begin load_negative(A[i,k]); multiply_by(c); add(A[j,k]);
	store(A[j,k]);
	end
@* The main program.
This is where the computation begins and ends.

@p begin @<Init...@>;
while not eof do new_vector;
@<Print the remaining null space vectors@>;
end.

@ When there is no more input, it's time to output the vectors that
generate the null space (if any).

@<Print the remaining...@>=
write_ln('The following vectors generate the null space:');
for i←1 to n do if null_space[i] then
	begin write('(');
	for j←1 to n do
		begin write_fraction(A[i,j]);
		if j<n then write(',')@+else write_ln(')');
		end;
	end
@* Index.